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A previously developed field-theoretic model [R.D. Coalson et al., J. Chem. Phys. 102, 4584 
(1995)] that treats core collisions and Coulomb interactions on the same footing is investigated in 
order to understand ion size elTects on the partition of neutral and charged particles at planar in- 
terfaces and the ionic selectivity of slit nanopores. We introduce a variational scheme that can go 
beyond the mean-field (MF) regime and couple in a consistent way pore modified core interactions, 
steric effects, electrostatic solvation and image-charge forces, and surface charge induced electro- 
static potential. Density profiles of neutral particles in contact with a neutral hard-wall, obtained 
from Monte-Carlo (MC) simulations are compared with the solution of mean-field and variational 
equations. A recently proposed random-phase approximation (RPA) method is tested as well. We 
show that in the dilute limit, the MF and the variational theories agree well with simulation results, 
in contrast to the RPA method. The partition of charged Yukawa particles at a neutral dielectric 
interface (e.g air-water or protein-water interface) is investigated. It is shown that as a result of the 
competition between core collisions that push the ions towards the surface, and repulsive solvation 
and image forces that exclude them from the interface, a concentration peak of finite size ions sets 
in close to the dielectric interface. This effect is amplified with increasing ion size and bulk con- 
centration. An integral expression for the surface tension that accounts for excluded volume effects 
is computed and the decrease of the surface tension with increasing ion size is illustrated. We also 
characterize the role played by the ion size on the ionic selectivity of neutral slit nanopores. We 
show that the complex interplay between electrostatic forces, excluded volume effects induced by 
core collisions and steric effects leads to an unexpected reversal in the ionic selectivity of the pore 
with varying pore size: while large pores exhibits a higher conductivity for large ions, narrow pores 
exclude large ions more efficiently than small ones. 

PACS numbers: 03.50.De,05.70.Np,87.16.D- 



I. INTRODUCTION 



Ion specific effects are relevant to an important number of biological, industrial and interfacial systems [T]. To 
mention a few examples, such effects are believed to play a crucial role in the absorption of large ions onto the water- 
air interface [lHl]j colloidal interactions in electrolyte solutions [5], charge reversal phenomenon [6 , protein stability [7], 
and the ionic selectivity of biological transmembrane channels that may discriminate between two different species of 
the same charge [SHTO]. 

One of the earliest and probably the best known experiments that brings to light ion specific effects is the measure- 
ment of the positive surface tension of inorganic salt solutions by Heydweiller [llj . The observed enhancement of the 
surface tension of electrolytes compared to that of pure water was explained by Wagner |12j in terms of the screened 
image interactions of ions that lead to their depletion at the water-air interface. Within this theoretical picture, 
Onsager and Samaras obtained their celebrated limiting law that relates the surface tension of electrolytes to the bulk 
electrolyte density [13 . This limiting law is indeed able to reproduce experimental data for dilute electrolytes |14) . 
but fails at large concentrations [15]. It is now well established that the discrepancy is due to the insensitivity of 
the Onsager-Samaras theory to ion specific effects, which are known to come into play at high densities. An earlier 
experimental observation of ion specific effects was made by Hofmeister who ranked salts according to their ability to 
precipitate egg white lysozyme [T5] . It is known that there exists some correlation between Hofmeister series and the 
surface tension of electrolytes. Although it is believed that Hofmeister series are also correlated with the size and the 
intrinsic polarizability of ions |17j , a universal theory able to explain the chemical characteristics that characterize 
the series is still lacking. Generalized electrolyte models taking into account at a MF level ion specific effects, such 
as ionic polarizability [ST", dielectric decrement [22] and ionic hydration [53] have been developed. Various theories 
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taking into account steric effects [THl HH] and ion polarizability [2D] have been also proposed in order to improve on 
the Onsager-Samaras theory. 

The finite ion size effects play also a key role in the case of electrolytes in contact with a strongly charged protein. 
Experimental evidence for such systems clearly shows ionic saturation effects close to the protein surface [23], which 
originate from close packing of ions. A modified Poisson-Boltzman theory that introduces steric effects by distributing 
the ions over a lattice have been proposed in Ref. [25]. A different field-theoretic approach that includes excluded 
volume effects was also introduced in Ref. Mean-field analysis of these theories was shown to successfully 

reproduce the saturation effect in question. However, these statistical models do not include a well-known propriety 
of molecular liquids with repulsive interactions, namely the wetting of the interface by the fluid. 

The structure of heterogeneous fluids has been a very active research domain over the past 50 years [27ll29] and 
early theories based on the solution of integral equations and mean-spherical approximations (MSA) were shown 
to reproduce with great accuracy the density profile of hard-core (HC) liquids in contact with an impenetrable 
wall [Sni ED ■ Recent works have focused on the structure of heterogeneous liquids interacting with short-range square 
well potentials [32] , Lennard- Jones [33] and HC- Yukawa potential mixtures [34] within different perturbative density 
fimctional methods (DFT). The first analytical solution of the Ornstein-Zernike (OZ) equation for a purely repulsive 
bulk Yukawa liquid was derived within MSA in Ref. 35 . Using a MF level DFT approach, it was recently shown in 
Ref. |36j that a binary heterogeneous fluid composed of repulsive Yukawa particles exhibits a wetting on the hard wall. 
Using a different method that consists in solving the OZ equation with the RPA closure, an analytical expression for 
the density profile of a repulsive Yukawa fluid was derived in Ref. j37j . In contrast to Ref. [3S] , the theory predicted 
a non-monotonic density shape, characterized by a depletion layer next to the wall, followed by a concentration peak 
and a decreasing particle concentration with increasing distance from the surface. 

The coupling between electrostatic interactions and core collisions has been considered beyond the MF limit within 
integral equation theories. Using the reference hypernetted chain approximation, the role of various mechanisms, 
including ion size effects on the charge inversion phenomenon was analyzed in refs. |38ll39j . The authors investigated 
also the effect of image charge and ionic dispersion interactions on the interplate pressure and ionic partitions in 
Ref. j40| . These integral equation methods developed within the particle density representation are efficient and 
provide good agreement with MC simulations [41| . However, they are also known to be computationally intensive 
and hard to use if one wishes to consider geometries other than planes. An alternative field theoretic model that 
couples electrostatic and core interactions was proposed in Ref. |42|. In this work, the void around an ion induced by 
the ion size was modeled with a repulsive Yukawa potential. By passing from the density to the field representation, 
the partition function of the system can be recast into a functional integral over two fields, namely a Yukawa field 
associated with core collisions and an electrostatic field induced by Coulomb interactions between charges. The 
authors considered the MF limit of this model in order to investigate the effect of ion size on the ionic partition next 
to two charged spherical colloids as well as on the electrostatic interaction energy between the colloids. In the present 
work, we reconsider this field theoretic model within three computation schemes, among which a variational approach 
inspired by a previous variational method that has been developed for ions with purely Coulombic interactions at 
charged dielectric interfaces [43] . The variational formalism of Ref. fl^ was also applied to the problem of ionic 
exclusion from cylindrical pores and yielded a new type of liquid-vapor phase transition that we proposed as the 
underlying mechanism behind the rapid switching of nanopore conductivity observed in experiments [441 145] . 

This article is organized as follows. The derivation of the fleld theoretic model for charged Yukawa particles is briefly 
explained in Section [iTJ Section |III| deals with a bulk Yukawa liquid. Using a general variational ansatz, it is shown 
that at the first order variational level, short-range core interactions within the homogeneous medium experience a 
further screening but Debye-Hiickel interactions remain unchanged. In Section. IV the physics of a heterogeneous 
Yukawa liquid is investigated. Various computation schemes are introduced and their validity for the density profile of 
neutral particles in contact with a hard wall is tested by comparison with MC simulations. The variational formalism 
is then applied to charged Yukawa particles in order to investigate core collision effects on the ionic distribution at a 
planar dielectric interface. An integral form for the surface tension taking into account the contribution from excluded 
volume effects is derived as well. In line with the Gibbs adsorption isotherm, core collisions that push the particles 
towards the interface lower the surface tension of the charged fluid. The variational equations are also solved in a 
slit pore in order to understand ion size effects on the ionic selectivity of slit nanopores. It is shown that due to 
a complex coupling between steric and excluded volume effects, image-charge forces and pore-modified screening of 
Coulomb and Yukawa interactions, large ions are more efficiently excluded from narrow pores than small ions while 
the penetration of large ions over small ions is favored into pores at nanometer scales. Limitations of the variational 
scheme, potential improvements and applications are discussed as well in the Conclusion. 
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II. MODEL 



We present in this section the derivation of the field-theoretic partition function for charged Yukawa particles. Since 
the derivation was already introduced in Ref. 221, we will exclusively report the general lines of the calculation. The 
canonical partition function of ions interacting with an electrostatic and a core potential reads 

^ =11 ^ J n dx,e--^({-^»--«(i-^» (1) 

1=1 T j = i 

where p is the number of ionic species, Ni is the number of ions for each species and At is the thermal wavelength of 
an ion. The electrostatic and repulsive Yukawa interactions are respectively given by 

Hc{{m,}) drdr'pc(r)«c(r,r')pc(r') (2) 

Hy{{^.,}) = \j drdr'pp(r)«;(r,r')pp(r') + j dvV^{v) p^^iv) (3) 

with the charge density 

P Ni 

Pc{'c)^^^(l^^i^-^^l) + '^s{r) (4) 
where qi stands for the valency of mobile ions, crs(r) is a fixed charge distribution, and the particle density is given by 

P Ni 
z=l i = l 

The Coulomb and the short-range Yukawa potentials are defined respectively as the inverse of the following operators 

v-\v,v') = -M^V[e(r)V<5(r-r')] (6) 
w~\vy) = ^^5(r-r') (7) 

where A = is the Laplacian operator. We introduced in Eq. ^ a spatially varying dielectric permittivity e(r). 
The Yukawa potential is obtained by inverting Eq. (|7| in Fourier space, which yields w[r) = ^j,e~''l''l/|r|. The wall 
potential Vw (r) in Eq. ([s]) takes into account the interaction between the particles and the boundaries of the system. 
It can be used either to model a specific particle-wall interaction, or to restrict the position of particles within a 
particular region. The self-energy of ions that has to be substracted from the total Hamiltonian is given by 

= |^;^^(r-r')|r=r' + ^z«(r-r')|r=r', (8) 
where the bare Coulomb potential v\{v — r') is the inverse of the Laplace operator 

-r'(r,r') = -^^A5(r-r'). (9) 

The electrostatic potential, solution of Eq. ([9|, has the well-known form wj(r) = £^/|r|, where = / [ATiewkBT) ~ 7 
A is the Bjerrum length in water solvent at temperature T = 300 K, e^, = 78 Eq stands for the dielectric permittivity of 
water and e denotes the elementary charge. The partition function of the system can be put in a more tractable form 
by performing two Hubbard-Stratanovitch transformations in order to linearize the pairwise interaction terms within 
the exponential of Eq. ([T]) at the cost of introducing two fluctuating auxiliary potentials for each interaction, namely 
an electrostatic potential </> and a Yukawa potential ?/;. The grand canonical partition function defined according to 
Zg — nr=i TliNiyQ e^^^'^Zc can be recast in this way in the form of a functional integral over these two potentials |42j, 

Zg= [ e-f^['^''^l (10) 
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with the functional Hamiltonian 

l2 



/ dr 



[V0(r)]^ 



87r^i3(r) 



ias{r)(j){r) 



dr 

8^ 



[V^(r)]V 62^2(^)1 _^A, / dre^^-^^W+'I^-^W+'^Wl. (11) 



In Eq. (11), we introduced the spatiaUy varying Bjerrum length ^^(r) = e / [47re(r)fcBr] and the rescaled particle 

4. 



fugacity Xi — e^' /X^ 



III. BULK IONIC LIQUID 



In this section, we analyze the electrolyte model of Eq. (Ill within the variational formalism in a bulk medium. 



The absence of boundaries in the system means a vanishing wall potential T4j(r) = 0. By taking advantage of the 
spherical symmetry within the homogeneous medium, we first introduce the Fourier-transformed potentials in the 
form 



Vc{r - r 



(27r)3 



>(g), ^(r) 



(27r)3 



iq-(r-r')~ 



(27r)3 



Vciq), w{r- 



w{q) 



with Vc = 4:TT£^/q^ and w = ATily/{q^ + 6^). By injecting these expansions into the total Hamiltonian Eq. (11) with 
^B(r) — (-w — const, one obtains 



H = 



1 



d^q 

(27r)3 



^(9)S-^(g)0(-9) + ^(<z)^-^(<z)^(-9)] -EV.e^^ / dre^^^[^-^(')+'^(^^]^"''-^ 



In order to simplify the field theoretic calculation that will follow, we re-express this Hamiltonian in a matrix form 



H 



1 



(27r)3^ 



dre*-' (2")' 



(12) 



where the q-dependence of the functions within the integrals in Fourier space is implicit. We have defined above the 
following matrix and vectors 



V = 



ic{q) 
w{q) 



# = ( til) 



Because the Hamiltonian functional (12) is non-linear in cj) and tp, an exact evaluation of the partition function is not 



possible. We will thus find the upper boundary to the exact grand canonical potential by making a variational ansatz. 
We choose a trial Hamiltonian of the most general quadratic form 



with the variational kernel 



d3q 
(27r)3 



G 1 = 



G~— 1 1 




G ^ G 



(13) 



(14) 



The reference Hamiltonian (13) contains a variational external potential vector Jv that will be defined below as 
well as the inverse of three trial potentials : the electrostatic potential G'o(r — r'), the core potential Gj,(r — r') 
and finally the operator G'c(r — r') that couples the fluctuating electrostatic and Yukawa potentials. The variational 
grand potential that will be optimized with respect to these trial functions reads Vli = JIq + {H — i?o)oi where 
f^o = - In/ e-'f^«["^''^l and the bracket (•)o means that the statistical average should be taken with the reference 

Hamiltonian Hq. A standard field-theoretic calculation yields 



/i = 



2 J (27r)3 ^"^^ ^ V J 



(15) 
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where we introduced the following functions 

F{q) = Tr (l - GV-i + In g) + J^GV-^GJ^ 
H,{q) = Jf GJ, + J^GJ, + Jf GJv. 



(16) 
(17) 



I stands for the 2x2 identity matrix and the trace operator acting on a general matrix M of the same size is defined 
as Tr(M) — J2n=i 2 ^nn- Let us note at his stage that the average value of the external fields computed with the 

Hamiltonian (13) is given by = (^(f)l,ijj^^ = (^^'^'^ — ^JvG. By injecting the inverse of this relation that gives 



Jv in terms of the variational potentials and Eq. (14) into Eqs. (16) and (17), one obtains 



Fiq) - 2 
H^iq) 



G— 1/^—1 
^ 

— 1 — 1 2 

^ 



In ( G^^G-i 



V 



■ipb{q) + qi4>b{q) 



(18) 



(19) 



The external electrostatic and Yukawa potentials follow from the variational equations dfi/S(t>biq) = and 
6 fi / Sipbiq) = 0. By taking into account the electroneutrality condition J2i Pb.iqi — 0, these equations yield 



E 



Pb,i 



(20) 
(21) 



where 



Pb,i = -A 



b.i 



dh 



b,i 



(22) 



is the average density of the charged Yukawa liquid. The remaining variational equations SJi/SGq ^ = 0, Sfi/SGy ^ = 
and 5fi/5G~^ = whose solution form with Eqs. (20) and (21) an upper boundary to the exact grand canonical 
potential now read 



y 

-2 



J + v, 



,-.-1 



Gn 



G ^ 



PtOt +W ^ +Gy^ 



Ptot + W ^ 



J + v7^+Gn^ 



G-i{g-2 + G-i [j + i-^]+Go' 



Ptot +W ^ - Gy ^ 



= 



with the total density ptot = X)i Pb,i ^^'^ the ionic strength J = Pb.iqf- With a little bit of algebra, one finds that 
the only solution to the above system of equations is 

G7^ = 



Gq ^ — 



g: 



q^ + «Dg 

q^ + 



yb 



AttL, 



(23) 
(24) 

(25) 



where we have introduced the Debye-Hiickel (DH) screening parameter ~ Airi^ ^ • pb.iqf and defined a new 

screening parameter 



(26) 



for the Yukawa potential. 

We note that by expanding Eq. ( [IT| ) at the quadratic order in (j) and ■0, one can see that at the level of the Gaussian 
theory, the electroneutrality condition Pb,iqi = imposes a vanishing coupling between these two potentials. 
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FIG. 1: (Color online) Geometry for a slit-like pore of thickness d. The dielectric permittivities of the pore and the substrate 
are and Em , respectively. Ions occupy the region < z < d — au,, where is the width of the Stern layer. 



According to Eq. (23), this coupling is also absent at the variational level. Furthermore, Eq. (24) indicates that 



the electrostatic potential in the bulk system is identical to the DH potential, that is, it is not modified by core 
interactions. But the HC potential given now by Eq. (25 ) gets an additional screening from the surrounding particles. 



It is important to note that in contrast vifith the interactions of electrostatic origin which are screened by the charge 
density, core interactions are screened by the particle density and this screening survives in the case of neutral particles 
{qi = 0). In the next section, we will show that the modification of this screening at the boundaries has important 
consequences on the distribution of particles close to solid interfaces. Finally, we note that by taking into account the 
equations ( 23p5 ), the relation Eq. (22) between the particle fugacity and the bulk density becomes 



Pb,i 



= A 



b,ie 



(27) 



IV. SLIT GEOMETRIES 



In this section, we investigate the effect of finite ion size on the configuration of ions at planar dielectric interfaces 
as well as on the ionic exclusion from slit pores. The slit pore of thickness d is characterized by a piecewise dielectric 
permittivity profile 

e{z) = e^e{z)e{d - z) + £™[0(z) + 9{z - d)], (28) 

where 9{z) stands for the Heaviside distribution and e„, and Sm are the dielectric permittivity of water and membrane, 
respectively (see Fig.jlJ. The pore is filled up with charged Yukawa particles whose confinement geometry is imposed 
with a steric potential defined as 

Vu,{z) — , ayj < z < d — aw 

Vw{z) — GO , z < Qw and z > d ~ a^, (29) 

where we introduced the width of the Stern layer Uw (the same for all ionic species) corresponding to a distance 
of minimum approach to the interface. In this work, we will consider exclusively the case of symmetric electrolytes 
composed of two species, with pI = Pb ^ Pb ^^i^ = —q^ = q, where + is for cations and - for anions. Furthermore, 
the pore is in contact with an external particle reservoir at the extremities and the particles within the pore are at 
chemical equilibrium with the reservoir, which fixes their fugacity according to = At^i. One thus obtains with 



Eq. (27) the equation that relates the fugacity of the particles within the pore to their reservoir density in the form 

Ai = pj, .e^'-^('^«''-'')-4^"''^^«. (30) 
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The parameters of the Yukawa potential can be indeed determined from neutron-scattering and X-ray experiments 
that measure the radial distribution function of liquids [29) . In this article the model parameters will be fixed from 
simple energetic arguments The dissipation length characterizing the range of the core interactions is usually 

chosen as a fraction of the ion diameter [ISl I47j. We will thus take b^^ — Ui, where stands for the ion radius. 
Furthermore, the strength of the Yukawa potential will be determined such that when two Yukawa particles are in 
contact, the interaction energy is on the order of the thermal energy, i.e. Vy (|r — r'| = 2ai) — A;bT/2, which fixes the 
intensity of the Yukawa potential according to £y — e^a^ ~ 7.39 a^. With this requirement, a closer approach of ions 
will be energetically unfavorable. 



A. Calculation schemes 



The field theoretic model Eq. ( 11 ) was studied in Ref. |42j for two charged spherical colloids immersed in a counterion 
liquid. The authors solved the MF equations with a lattice MC procedure. Correlation effects neglected at the MF 
level can be partially captured with the RPA closure of Ornstein-Zernike (OZ) equations, introduced in Ref. [37j . 
Below, we will first explain the derivation of the MF equations and the associated boundary conditions for the general 
case of charged particles and boundaries carrying a fixed charge distribution. We will then report the prediction 
from the recently developed RPA theory for the same system [37]. An alternative calculation method consists in 
making a variational Ansatz for the reference Hamiltonian and finding the optimal choice for the trial parameters 
from the minimization procedure. Density profiles obtained within these calculation schemes will be compared with 
MC simulations in Sec. |IVB| for the case of neutral particles in contact with a planar wall. 



Mean-Fteld limit and RPA approach 



The mean-field (MF) solution of the grand canonical partition function Zq associated with the weak surface charge 
and low density limit is obtained from the saddle point equations SH/6(f){r) = and 5H/5%1}{t) = [H], which read 



A(/.(r) - 



g-y„(r)-^(r)+^, sinh(/)(r) = -47rqCcr,(r) 
&V(r) + 87r^yPbe-^-('-)-^'("-)+^"' cosh0(r) = 0, 



(31) 
(32) 



where 



"DH 



j^Pb, ij'b — Sniy/b'^ and we have also defined ip{r) — —iip{r), 0(r) — —iq(f>{r) and dropped the bar 



over the potentials in order to simplify the notation. In the case of a charged HC liquid confined between two planar 
surfaces located at z = and z = d, with surface charge a{z) — CTs [S{z) + S{z — d)], the differential equations (ISll) 
and ( 32 1 should be solved in the region a^u < z < d — with the following boundary conditions (see Appendix [ 



dz 

d^ 
dz 

d4 
dz 

di! 
dz 



btp {z 



z — {d—aT^ ) 



z — {d~au; ) 



-bip (^z = (d- a^) ) 



(33) 
(34) 
(35) 
(36) 



In the case of a system composed of a single interface separating a particle free half space and a solvent medium filled 
up with ions, the boundary conditions that replace Eq. (35 1 and (361 read 



(j){z oo) = 



(37) 
(38) 



where the constant bulk value of the Yukawa potential is given by Eq. (21) 



In this article, we will treat exclusively the case of neutral interfaces, i.e. ct^ = 0. In this limit, the external 
electrostatic potential 4'{z) vanishes everywhere. In the dilute limit ^irlyPb < b^, by expanding at the linear level in 
^{z) — ipb the remaining MF equation (32), one can easily find a linear solution for the Yukawa potential ^{z). The 
equation to be solved reads 



1^ 



(39) 
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The solution to this equation satisfying the boundary conditions ( |34[ ) and (|36| is given by 



hijji, cosh [Kyb [d/2 — z)] 



Kyb sinh [Kyii{d/2 — a^)] + 6 cosh [Kyi,{d/2 — a^)] 
For the single interface system that can be recovered in the limit d — > oo, this linear solution reduces to 

b 



(40) 



-Kyh{z — a^) 



^yb 



(41) 



Finally, the density profiles are obtained by injecting the potentials (40) or (41) into the MF relation 

p{z) ^ pte^^-^^'\ 

In a recent work, particle density profiles were derived from the solution of the OZ equation with RPA closure 
for the single wall system. The RPA result reads 

Prpa{z) = pb[l- Vy{z) - 5i^RPA{z)] 

where the external potential reads 

5^rpa(z) = -tpRPAiz) - A-KlyPb/b^ = -4TT£yPb/b'^e~'"'. 



(42) 



(43) 



(44) 



We note that the external potential ipRPAiz) automatically satisfies the boundary condition Eq. (34) on the wall. 
Although the RPA method partially accounts for correlation effects absent at the MF level, we will show in Sec. |IVB| 
that the MF theory provides a better agreement with MC simulations. 



2. Variational approach 



An alternative approximative way that allows to go beyond the MF level consists in making a variational ansatz 
in order to compute an upper boundary to the exact grand potential Q = —ksTlnZQ. In Sec. Ill where we treated 
the case of a bulk electrolyte, it was shown that at the variational level, the electroneutrality condition imposes a 
vanishing coupling between the electrostatic and Yukawa potentials. Although this result does not imply a vanishing 



coupling within the pore, a variational Hamiltonian of the form (13) that couples 4'{r) and ijj{r) would lead to a 
very complicated and untractable set of variational equations for the pore system. We thus opt for a restricted trial 
Hamiltonian that treats (f>{r) and ip{r) separately, namely Hq = Hq^ + i/o0 with 

Ho<i,^l I [m-^M^)]vo\ry)W)-^M^')] 



r.r 



H^i' = \ / [^(r) - *^o(r)] w^\v, r') [^^(r') - #o(r')] 



r.r' 



(45) 



The reference Hamiltonian (45) contains four trial functions : two external potentials (/'o(r) and "00 (r) whose physical 
origin will be clarified below, an electrostatic potential vq{y,y') and a core potential Wq^t^y'). The variational grand 
potential reads fiy = fio + {H — Ho)^, where fig = —IuZq = — In/ Vcj/Dt/) e~^°^'l'''^^ is the gaussian part of the 
grand potential that corresponds to the vdW theory [49]. Since we did not consider an explicit coupling in the 
reference Hamiltonian Hq [0, ip] , the gaussian contribution can be separated into the Coulomb and Yukawa parts as 
flo = 51o0 + f^oV'i where 

f7o0 = - In Jv(t> e-^°*["^l (46) 

Oov = ~ e"^"*!'^!. (47) 

Evaluating the functional integrals in Qy, the variational grand potential becomes 

= no4, + no^+ J dr Ul^^^ +a3(r)0o(r) j - J ^ ([V0o(r)] V ^Vo^W) 



J drdr'<5(r-r')Vr£(r)Vr'Z;o(r,r')+ J ^(5(r - r') (VrVr^ + 6^) wo(r, r') - J] J drp^r) (48) 
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where we defined the local density as 



(49) 



The variational equations for the trial potentials i5ri„/(5(/)o(r) = 0, Snv/6ipo{r) — 0, 6^v/6vo{r,r') — and 
S^lv/Swo{r,r') = read 



A0o(z) + 47r£s(r)2j 

Pi{r)qi = -47r£B(r)(Ts(r) 







(50) 
(51) 



[-V(£(r)V) + e{v)Kl{r)] vo{r, r') = — 
[-A + K.lir)] wo{r, r') = Aniydir - r') 
where we introduced the spatially varying screening parameters as 

'^c(r) 



<5(r- 



47r^y^p,(r). 



(52) 
(53) 

(54) 
(55) 



The coupled self-consistent equations (50)- (53 1 take into account correlation effects at a non-linear level. Eqs. (50) 
and ( 52 ) were derived within a variational calculation without Yukawa interactions in Ref. [48^ . Eq. ( 50 ) is an extended 
PB equation whose solution yields the electrostatic potential in the slit pore, dressed by pore-modified correlation 
effects. Eq. (51 ) provides the local value of the external Yukawa potential tpoiz) that quantifies excluded volume effects 
in the pore. The third relation Eq. (52) is a modified DH equation and Eq. (53 1 takes into account the modification 
of Yukawa interactions in the slit pore. These type of non-local closure relations are known to be unsolvable even in 
the restricted case of ions without core interactions [121 135] . Hence, to make further progress we restrict the form of 
the pairwise potentials vo{r,r') and wo(i",i"') as follows. 

Following our work on simple ions [33], we choose the variational electrostatic potential as the inverse of the 
generalized DH operator 



^(r,r') 



[-V(e(r)V) + e(r)K2(r)] 5{r - r') 



(56) 



where £(r) is given by Eq. (28) and Kc(i") — Kc9{z — aii;)9{d^ 



z) is a constant piecewise trial screening parameter 



which is uniform within the pore. This choice is indeed a generalized Onsager-Samaras approximation [13J. It was 
shown in our previous work on purely Coulombic liquids that the approximation yields a good agreement with MC 
simulations for ion distributions in charged pores 53] • 

For the sake of physical consistency, the choice of the trial Yukawa potential should be done with care. In section HI, 
it was shown that Yukawa interactions get an additional screening due to the presence of surrounding particles. Hence, 
we expect the inverse screening length to reduce to the bare value b in the particle-free regions. This point is also 
confirmed by the form of the screening parameter Eq. (55) associated with the general variational equations (50)- (53 1. 
We thus chose the variational Yukawa potential as the inverse of the following operator 

-A + nlir) 



Wn 



-S{r- 



(57) 



where Ky{r) = b[9{aii, — z) + 6(z + a^u — d)] + Ky9(z — aw)d{d — — z) is an effective screening parameter characterizing 
the pore-modified screening of core interactions that reduces to the bare screening parameter in the membrane. 

The effective electrostatic and Yukawa potentials obtained by inverting Eqs. (56) and (57) in the slit geometry are 
reported in Appendix|Bj We are thus left with two variational parameters Kc, Ky and two variational potentials (f>o{v) 
and V'o(r). The corresponding variational grand potential reads 




b'roiz) 



dzp,{z) 



(58) 



(59) 
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where S stands for the lateral surface. For charged Yukawa particles confined in a slit pore, by taking into account 
the chemical equilibrium condition Eq. ( 27 1 and defining the following potentials 



^c{z) = y [iw{KDH - lie) + 5vo{z)] 



the local particle density Eq. ( 49 ) becomes 



(60) 
(61) 

(62) 



For a symmetric electrolyte, by rescaling the electrostatic potential according to 4'{z) = q<j){z), the variational equations 
for the external potentials, S^ly / 5iI]q{y) — and 5r2t,/(5^o(r) = read 



A(/.o(z) 



-K„(^)-i/o(^)-i/,(^)+V6-Vo(2)sinh0o(^) = "47rC'?^Ts(2) 



A^o(-2) - ''Vo + 87r£j,pf,e-^™(^)-^'^(^)-^''(^»+'''''-'^°W cosh<?!>o(z) = 



(63) 
(64) 



where we dropped the bar over the electrostatic potential in order to simplify the notation. The differential equa- 

The variational equations for the re- 
and dflv/dKy 



tions (63) and (64 1 should be solved with the boundary conditions (33)-(38) 



maining trial parameters can be obtained in a similar way from the simple derivatives dfly/dK^ 
in the form 

\ / p \ ^y I p 

where we have introduced the pore average defined as = dz ■ / {d — 2aw). 





(65) 
(66) 



The set of coupled variational equations ( 63 )-( 66 1 is the main theoretical result of this article. The relation Eq. ( 65 ) 
is a modified DH equation for the effective pore screening length. Eq. (66) takes into account the modification of the 



screening of Yukawa interactions in the slit pore. It is important to note at this stage that these relations are valid in 
all geometries (spherical, cylindrical ...) and can be solved with an iterative algorithm. We emphasize that Eqs. (63) 
and (65) were derived for ionic liquids without Yukawa interactions {tp{z) = and Vy{z) = 0) and solved in slit and 
cylindrical pores in Refs. [43 45. . It should be also noted that approximate self-consistent equations, similar in form 
to Eqs. (63) and (65) have been frequently used in nanofiltration studies [SUHS^] . The prediction of these approximate 



equations for salt rejection from neutral and charged pores was compared with that of the variational equations in 
Ref. [43[ 145) and it was shown that the former equations can overestimate ionic concentrations in cylindrical pores by 
a factor of 2 to 3. The relations Eqs. (63)-(66) now generalize these self-consistent equations to the case of electrolytes 



composed of large ions associated with important excluded volume effects. 

The potentials (60) and (61) are obtained in the limit r' — >■ r from the corresponding kernels given by Eqs. (BIO) 



and (B16). One finds for the slit pore geometry 



Veiz) 



-^{KdH - Kc) + - 



2 

dkk 



dkk 



Pc 1 - A2e-2Pc(rf-2a„) 

A 



-2pc(2-a„) _|_ g-2pc(d-a„-z) _j_ 2^^g-2pc((i-2a„) 



g-2p„(2-a„) _|_ g-2p„(d-a„-2) _|_ 2/\^g-2pa (d-2a„) 



(67) 
(68) 



The functions pc, Py,ric, Ac, Ay are reported in Appendix [B} We also introduce the partition coefficient, i.e. the 
pore-averaged particle density as fc = (p(z)/pb)p. 

In the case of a simple dielectric interface located at z = 0, the minimization of fl^ at constant fugacity with respect 
to the inverse screening lengths yields Kc = hdh and Ky — Kyi,. We are thus left with two variational equations (63) 
and (64) to be solved. For the single interface geometry {d — > oo), the Yukawa potential (68) becomes 



dkk 
Pv 



(69) 
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and the electrostatic potential (67) that we separate into a solvation and a dielectric part, Vc{z) = Vcs{z) + Vci(z) 
reads 



2 J pc Pc + k 

J Pc {Pc + k){pc + rick) 



(70) 
(71) 



We note that in the case of neutral interfaces (cts = 0), the external electrostatic potential vanishes, i.e. (j)o{z) = 0. 

We also computed the surface tension for the charged Yukawa fluid. The surface tension, denoted by (Te, is defined 
as the excess Grand potential, i.e. the difference between the Grand potential of the interfacial system and that of 
the bulk system. tJe was computed in Ref. W within a first order perturbation theory. The model included a Stern 
layer induced by the finite ion size. A sophisticated variational calculation of CTe that accounts for a depletion layer 
induced by electrostatic image forces was also proposed in Ref. The derivation of from the variational Grand 
potential Eq. ( 58 1 is explained in Appendix [C] The result reads 



247r 



yb 



■° dz 



^DH\ 



dfcfc 
"8^ 



In 



{pc + kricY kI,jjA 



8nL 



[(VV'o)' + 6'(V'o-^?)] 



kpciVc + 1)2 

DO 

dz [p{z) ~ pt] 



2p? 



(72) 



where /? = l/fc^T. The relation ([72]) generalizes the result of Ref. [19 to the case where core collisions are present. 
The first three terms proportional to include the contribution of the exclusion effect induced by the Stern layer. 
The next four terms contain corrections from the first order variational calculation to the quadratic fiuctuations of the 
electrostatic and Yukawa potentials. Finally, the last term brings non-linear contributions from the surface depletion. 
We finally note that the evaluation of Eq. ( 72 ) requires the numerical integration of Eq. ( 64 ) . 



B. Single interface 



In this section, we investigate the partition of neutral and charged Yukawa particles in contact with a neutral planar 
interface located at z = 0. The single interface system is recovered from the slit geometry depicted in Fig. [l] in the 
limit d — > (X). The right half space is filled up with ions in water medium of dielectric permittivity = 78. The 
left half space is free of ions and composed of biological matter or air, associated with a low dielectric permittivity 
£m < £w We will take Em = 2 (the case for lipid membranes), unless otherwise stated. This choice is also a good 
approximation for the water- air interface characterized by Sm = 1- We will first investigate the case with a vanishing 
Stern layer a^, — and then discuss the modifications induced by a finite Oy, . 

Fig. 2 displays for Uw — the distribution of neutral Yukawa particles {q — 0) for model parameters Xy = = 

1/3 

2.992 and Ty = (.yPf^ = 0.25, 0.4, and 0.6. We present a comparison of the MF theory, the variational formalism and 
the RPA method with the result of canonical MC simulations. First of all, one sees that MC data show wetting of the 
interface. The adsorption of repulsive particles onto a hard wall was also observed in other MC simulations [5^IM1[51] 
as well as within various theoretical approaches [311 1371 155j . The underlying mechanism is known to be the pressure 
applied by the bulk on the particles close to the surface, which originates from core collisions. At the MF and 
variational levels, this effect is incorporated in the attractive potential ~ 4'b obtained from the numerical solution 
of Eqs. ( 32 ) and ( 64 ) . By comparing the plots of Fig. [2] from top to bottom, one sees that an increase of the coupling 
of core interactions that amplifies excluded volume effects results in a stronger particle adsorption on the wall. A 
comparison of the MF level density profile with the simulation data in Fig. [2] shows that although there exists a good 
agreement far from the interface, the MF theory slightly overestimates the particle density close to the surface and 
the discrepancy grows with increasing Ty. We also note in passing that the analytical solution given in Eq. (41) of 
the linearized MF equation in Fig. [2ja shows a very good agreement with the numerical solution of Eq. (32). 

The overestimation of the particle density at the MF level can be explained by the incapacity of the MF theory 
to take into account the correlations associated with core interactions that are partially embodied at the variational 
level in the potential Vy{z) of Eq. (69 1. In Sec. Ill, it was shown that in a homogeneous medium, a Yukawa particle 
surrounded by other particles experiences core interactions reduced by the screening effect. The repulsive potential 
Vy (z) that reduces the contact density is induced by the modification of this screening close to the boundary : because 
the left side of the interface is free of particles, the Yukawa particle is more efficiently screened in the bulk than in the 
neighborhood of the surface. Since the screening is favorable to the system, the particle exhibits a tendency to move 
away from the interface. As can be seen in Fig. [2j the variational formalism that can account for this weak solvation 
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effect shows a better agreement with MC results. One notes in Fig. [2]b and c that very close to the interface, the 
variational result deviates from the simulation data and exhibits a weak concentration peak absent in MC simulations. 
This might be due to the simple choice of a uniform Ky in Eq. ([57]). We emphasize that for a bulk concentration 
= 1 M, the region where the disagreement takes place corresponds to the interval z < 0.5 A. Finally, we note that 
the RPA method exhibits a poor agreement with MC simulations over the whole interfacial area and the discrepancy 
increases with Ty. 

We illustrate in Fig. [sjthe density of monovalent ions (q — 1) in the neighborhood of a dielectric interface (e.g. air- 
water or protein-water interface) for two bulk concentrations and four ion sizes. We note that in the case of particles 
with a finite charge, Eq. (64) is integrated by including the electrostatic potential Vc{z) given by Eq. ([70| and (71 1. 



This potential contains image-charge interactions induced by the dielectric discontinuity through the interface and 
an electrostatic solvation effect. The latter originates from the absence of ions in the region z < that modifies 
the electrostatic screening close to the boundary. Both effects are known to be repulsive and in the case of point 
like ions without core interactions, they lead to an ionic exclusion layer at the interface and a density profile that 
monotonically increases towards pb with increasing distance from the surface [431 153] . This behaviour is illustrated in 
Fig. [sjby solid black curves (a^ = 0). It is well established that the same image forces contain the leading contribution 
to the positive surface tension of electrolytes at the water-air interface [13j and the exclusion of ions from confined 
pores [331 mi inni IS3] • By comparing the density profile of point-like (a^ = 0) and finite size ions (a; > 0), one notices 
that core collisions lead to a net increase of the ion density in the interfacial region. For ions of size a.; < 2 A, the 
density profile p{z) keeps its monotonic trend whereas for larger ions (a^ > 2 A), p{z) exhibits an oscillatory shape. 
With increasing distance from the surface, p{z) exceeds the bulk density and exhibits a maximum at a characteristic 
distance. Beyond this distance, the ion concentration begins to decrease towards the bulk limit. Moreover, we 
notice in Fig. [3] that ionic adsorption becomes more pronounced with increasing bulk density. Specifically, for a bulk 
concentration pi, = 1 M, the amplification of the concentration peak that accompanies the increase of the ion size 
from = A to = 3.0 A is four times higher than in the more dilute case Pb = 0.3 M. 

As illustrated in Fig. |4] that shows the behaviour of the electrostatic and Yukawa potentials as well as the evolution 
of the total potential of mean force (PMF), the non- monotonic shape of the ion density originates from a competition 
between the combined image-charge and solvation forces (included in the potentials Vc{y) and Vy{z)) that drive the 
charges to the bulk and the positive pressure applied by the bulk particles on the surface particles (the external 
Yukawa field tpoiz)) that move them towards the surface. A comparison of the electrostatic and Yukawa potentials 
for the cases = 1.5 A and Ui — 3.0 A shows that the intensification of the ionic adsorption with increasing ion size 
is essentially due to a rise in the amplitude of the attractive potential ipiz) — V'fc, which in turn originates from an 
amplification of excluded- volume effects. Although the repulsive potential Vy{z) increases with increasing almost by 
the same amount as V'(^) ~ipb on the surface, we notice that it experiences a stronger screening than the latter. Hence, 
the effect of the potential Vy{z) comes into play over a short distance. Furthermore, by comparing the amplitude of 
the potentials in Fig. |4j one notices that in the limit = that we have considered so far, the ionic depletion is 
mainly due to the interaction of ions with their images, i.e. Vc^i{z). This is also illustrated in Fig.[3](b) where we plot 
for two particle sizes the density profiles for the case Em = £w = 78, where image interactions vanish (dashed lines). 
It is seen that the depletion layer is greatly reduced and the height of the concentration peak for the ion size Ui — 3 
A is amplified. 

We illustrate in Fig. [sja the evolution of the surface tension Ce Eq. (72 1 as a function of the bulk concentration 
for various ion sizes and the model parameters Sm = 2 and aw — 0. According to the isothermal Gibbs adsorption 
equation dag = —jdp,, where p is the chemical potential of the particle and 7 = dz [p{z) — pb] denotes the surface 
excess, the surface tension of a repulsive interface grows with pb whereas for an attractive surface, is a decreasing 
function of pb- The former case, associated with the primitive model of ions at the water-air interface |12| is illustrated 
m Fig.[5|a. The ionic depletion induced mainly by image-charge repulsion leads to an almost linearly increasing surface 
tension with pb. The first important point to note in this figure is the reduction of de with ion size. This feature is 
particularly noticeable for sizes = 2.5 and 3 A. The decrease of with increasing is due to the enhancement of 
the surface wetting of ions induced by core collisions (see also Fig. [3]) . Moreover, one notices that although Ue keeps 
its monotonic shape for = 2.5 A within the considered submolar concentration range, the curve for the largest ion 
size Oi = 3 A exhibits a maximum at ~ 0.7 M and decreases beyond this value. This reversal is also explained by 
the enhancement of excluded volume effects at large enough bulk concentrations and ion sizes, as confirmed by the 
comparison of Figs. [3] a and b. 

We also computed ion densities and the corresponding PMFs in the presence of a Stern layer with — ai. 
The result is displayed in Fig. [6)a. A comparison of Figs. |3] and Fig. |6] shows that besides a sharp cut-off of ion 
concentrations at z = a^, the most significant effect of the Stern layer is a decrease in the intensity of the repulsive 
image potential, which in turn leads to higher concentration peaks. This reduction of the image force is explained by 
the fact that the larger a^, is, the larger is the separation between the central charge and the interface. Indeed, one 
can check in Eq. (71 1 that the amplitude of Vci{z) exponentially decreases with increasing a^- A comparison of the 
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surface tension curves in Figs.jSja and b confirms this observation : for a given ion size a^, image forces are weakened 
by a finite a^), which results in a stronger particle concentration at the interface and lowers the surface tension. At 
this point, we emphasize that our mapping between the ion size and the parameters of the repulsive Yukawa potential 
is not unique. A different choice will modify the range of the core interactions and lead to a quantitative change of 
the curves in Figs, [sja and b. However, this does not change our main conclusions here. 

It was shown in this section that the consideration of core interactions between particles gives rise to a significant 
ionic adsorption at the neutral dielectric interface and the adsorption effect is amplified with the increase of the ion 
size or the bulk concentration. We also found that in agreement with Gibbs adsorption isotherm, this partial wetting 
lowers the surface tension of the primitive electrolyte model. We investigate in the next section the effect of core 
collisions on the exclusion of neutral and charged particles from slit pores. 



C. Slit pores 



In this section, we investigate the role of particle size on the exclusion of Yukawa particles from neutral slit pores, 
in contact with an external particle reservoir at the extremities. The pore geometry is depicted in Fig. [l] In the 



neutral particle limit q — where Eq. (651 yields Kc — 0, partition coefficients and density profiles are obtained from 
the iterative solution of the variational equations ( |64[ ) and (66) whereas for charged particles, the third variational 
equation ( |65[ ) for should be included as well in the iterative algorithm. 

Fig. [7|a displays the partition coefficient of neutral Yukawa particles of size at — 2 A versus the pore size for 
ttw = and two different bulk concentrations. One notices that the MF theory predicts a pore density higher than the 
reservoir density. The particle adsorption into the pore is explained by core collisions that lead to an accumulation of 
particles at the boundaries that prevent their expansion. Since the confinement becomes stronger with decreasing pore 
size, the pore partition coefficient increases in a monotonic way. While the prediction of the variational calculation 
remains qualitatively similar for large pore sizes, the evolution of k at small pore sizes significantly differs from the 
MF picture. At a characteristic pore thickness, k reaches a maximum and then starts decreasing. 

The reversal of the MF behavior originates from a competition between hard-core collisions and the pore-modified 
screening of these interactions. At the variational level, the latter effect is embodied in the integral term of the 



potential Eq. ( |68[ ) that we will denote by SVy{z). In Section IVB on single interfaces, it was shown that SVy{z) 
incorporates a solvation effect that pushes the particles towards the bulk region where Yukawa interactions can be 
screened more efficiently than next to the interface. This potential contains a similar effect in the slit pore. Indeed, 
dVy{z) excludes the particles from the pore. This feature can be understood by noting that particle penetration 
into the pore necessarily adds to the solvation and amplifies the intensity of SVy{z) (except at high densities where 
the screening of SVy{z) dominates the latter effect and SVy{z) begins to decrease with increasing ph). Moreover, the 
reduction of the pore size also increases SVy{z) since the modification of the bulk screening is more significant in small 
pores. Consequently, 6Vy{z) becomes more repulsive with decreasing d and at a characteristic pore size, its variation 
with d dominates that of the attractive Yukawa potential tpb — ipo{z) and k begins to decrease. A comparison of 
the curves for pi, = 0.2 M and pt, = 0.5 shows that an increase of the reservoir density that positively adds to the 
amplitude of these opposing forces leads to a stronger competition between them and a higher peak for k. We also 
show in Fig. [7]b the evolution of the variational inverse screening length Ky with the pore size. The inspection of 
the plot shows that Ky exhibits a similar shape to the partition coefficient, which is explained by the fact that Ky is 
an increasing function of the particle density in the pore (see Eq. (66 1). Finally, we note that in the case of neutral 
particles, adding a Stern layer simply reduces the accessible volume by 2Saw In other words, considering a finite 
Ou, is exactly equivalent to redefining the pore size according to d' — d ~ 2aw , which has the effect of shifting the 
curves in Fig. [7] towards larger pore sizes. This equivalence is valid for charged Yukawa particles exclusively in the 

We compare in Fig. [8] the density profile of neutral and charged Yukawa particles for aw = and two different pore 
sizes, namely d = 8 and 16 A. First of all, one notices that the pore density profile of neutral Yukawa particles exhibits 
an oscillatory shape, with two concentration peaks that originate from the competition between the pore-modified 
screening of core interactions and core collisions that push the particles towards the pore walls. While the density 
profile of charged Yukawa particles is also oscillatory for = 78, strong image forces that come into play in the case 



Em = 2 smooth the density profile. Furthermore, as in the case of a single interface that was investigated in Sec. IVB 
the density profile of finite size ions is monotonic, and reaches its maximum value in the middle of the pore [431 153j. 
This is explained by the fact that point-like ions are exclusively subject to repulsive image and electrostatic solvation 
forces, whose intensity reach their minimum value at z = d/2. More importantly, it is seen that even in the case of a 
vanishing dielectric discontinuity Sm = £w, the local density of monovalent ions is lower than that of neutral particles. 
The reduction of the particle density is due to repulsive electrostatic solvation forces embodied in the potential Vc{z) 
of Eq. (67) that acts exclusively on particles of finite charge. Moreover, comparison of the curves for = 2 in Fig.jsja 
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and Fig. [sjb shows that if the pore size is increased from d = 8 A to d = 16 A, the reduction of the intensity of 
image interactions leads to an extra ionic penetration into the pore. The underlying dielectric repulsion mechanism 
that is responsible for ionic exclusion from membrane nanopores was thoroughly investigated for point like ions in 
refs. [131[5ni[53]. 

We observe in Fig. [Sja and b that in the limit £,„ = £w, due to core collisions that push the particles into the pore, 
finite size ions have a larger density than point like ions. However, one also notices that the decrease of the membrane 
dielectric permittivity is accompanied with a highly interesting reversal. For e„ — 2, despite the core collisions that 
guide the finite size ions towards the pore, one finds that the local density of point like ions exceeds that of finite 
size ions. The same reversal effect is also illustrated in Fig. |9]a which displays the ionic partition coefiicient against 
the pore size for = 0.5 M and various values of Em- It is shown that in the case of a weak dielectric discontinuity 
and large pore sizes, the penetration of finite size ions is favored over that of point like ions, but for low values of 
that correspond to strong image forces, the former experience a significantly stronger exclusion from the pore. This 
seemingly counter-intuitive effect is mainly due to the extra solvation energy barrier for the penetration of finite size 
ions into the pore. The energetic barrier is associated with the screening of core interactions, whose contribution to 
the particle density corresponds to the term £y{Kyi, — Ky) of Eq. (68) that vanishes for Oj = (or £y =0). Indeed, 
a decrease of d or Sm, accompanied with an amplification of image-charge and solvation forces results in a reduction 
of the ion density and Ky (see Fig. [7]b). As a result, when the ionic exclusion becomes strong enough, which occurs 
for small pores or low values of em, the term ty{nyb — Hy ) reaches a large enough value to dominate the contribution 
from core collisions, i.e. the term ^{z) — ipb in Eq. (62) and the density of finite size ions drop below that of point 
like ions. 

Fig. [9]b displays the partition coefiicient of finite size and point like ions against their bulk concentration. It is 
shown that the reversal of the facilitated penetration of large size particles with decreasing e„ extends over the whole 
concentration regime, i.e. this peculiarity also survives in the very dilute limit. However, since the decrease of the bulk 
concentration results in a reduction of excluded volume effects, the difference between the pore density of point-like 
and finite size ions monotonically decreases with their bulk density. 

The evolution of variational inverse screening lengths for monovalent ions is illustrated as a function of the pore 
size in Fig.[7jb and Fig. [TOj First of all, we notice in Fig. [TO] that Kc also exhibits the reversal effect discussed above. 
Then, the inspection of Fig.[7]b and Fig. 10 shows that as expected from the variational equations (65) and (66), Kc 
and Ky qualitatively follow the same trend as the partition coefficient. Namely, in the case — where the ionic 
exclusion is not too strong, k^. and Ky monotonically decrease from their bulk values towards slightly lower values with 
decreasing d. In the opposite limit = 2 that corresponds to strong image forces, at a characteristic pore thickness 
below which the salt rejection is almost total ((p(z)/pb) ~ for d ~ 4 A), the system reaches the dilute limit, and 
one gets Kj, ~ & and Kj. — 0. 

We also plot in Fig. [TT] ionic partition coefficients as a function of the pore size with a finite Stern layer a„ = a-i- 
In addition to an effective reduction of the pore volume accessible to ions, it was shown in Sec. |IVB| that a finite 
increases the distance between the interface and the central charge, which results in an overall reduction of the 
image repulsion. Thus, one expects these two opposing effects to shift the curves in Fig.|9]a towards larger pore sizes 
along the abscissa axis and also to larger densities along the ordinate axis. As a result. Fig. [Tl] shows that for any 
membrane dielectric permittivity satisfying Sm < £io, ionic penetration into the pore is characterized by a stronger 
rejection of finite size ions at small d and their favored penetration over point like ions for large d. 



V. CONCLUSION 



In this article, we have studied a field-theoretic model that considers electrostatic and core interactions on an equal 
footing in order to investigate the role played by excluded volume effects on the partition of ions at dielectric interfaces 
and ionic selectivity of slit nanopores. To this end, we have developed a variational calculation scheme that allows to 
go beyond the MF regime and couples in a consistent way pore modified core interactions, steric effects, electrostatic 
solvation and image-charge forces, and electrostatic potential induced by a surface charge distribution. 

In the first part of the work, we considered a bulk liquid composed of charged Yukawa particles. Using a general 
variational ansatz, it was shown that at the first order variational level, short range core interactions between two 
test ions experience a further screening from the surrounding particles but electrostatic interactions between them 
are not modified by core collisions. 

The second part of the article was devoted to the configuration of ions at simple interfaces. Density profiles for 
neutral Yukawa particles obtained from the MF limit of the theory and the variational formalism, and also from a 
recent RPA method were compared with MC simulations. We showed that in the submolar concentration regime, the 
MF theory presents a reasonable agreement with simulation results, but slightly overestimates particle concentrations 
in the proximity of the hard wall. It was shown that the variational theory, which includes the correlation effects 
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associated with the modification of the screening of core interactions by the interface gives a better agreement with 
simulation results. We also observed a deviation of the variational prediction from MC results very close to the 
interface, namely for z < 0.5 A. This discrepancy is probably be due to the simple choice of a uniform variational 
screening length for core interactions that does not account for the variation of the particle density close to the wall. 
It was also shown that the RPA method overestimates concentration profiles over the whole interfacial area. In the 
case of charged particles in contact with a dielectric interface, the competition between core collisions and repulsive 
dielectric forces yields concentration peaks at the interface and the corresponding ionic adsorption is amplified with 
increasing ion size or bulk electrolyte concentration. We computed an integral expression for the surface tension that 
accounts for excluded volume effects and also a Stern layer, induced by the finite ion size. In agreement with the 
isothermal Gibbs adsorption equation, we showed that an increase of the ion size that yields a stronger wetting on 
the wall is also accompanied with a reduction of the surface tension. 

In the third part, we solved the variational equations in a neutral slit pore. We first treated the case of neutral 
particles to show that excluded volume effects associated with core collisions within the bulk reservoir result in a net 
particle adsorption into the pore. In the absence of a dielectric discontinuity, we found that due to the same core 
collisions that drive the particles into the pore, the density of ions increases with their size. However, the presence of a 
dielectric discontinuity reverse this picture : it was shown that the ionic selectivity of pores characterized with strong 
image forces increases with the ion size. The reversal of the pore selectivity originates from a complex coupling with 
the solvation energy associated with core interactions and image forces, and survives in the limit of dilute electrolytes. 
This interesting peculiarity calls for experimental confirmation. The present pore model includes as well a Stern 
layer associated with the surface affinity of ions. Besides a reduction of the volume accessible to ions, a finite 
decreases the intensity of repulsive image forces. As a result, small pores of subnanometric size exhibit a stronger 
selectivity for large ions but nanometric pores exclude small ions more efficiently than large ions. 

In this article, we have investigated exclusively the case of neutral planar interfaces. However, the variational 
equations (63)-(66) are derived for the general case of charged interfaces. The role played by excluded volume effects 
on ionic exclusion from charged pores will be presented in future work. We note that the relations Eqs. (63l-(66) 
can also treat the case of curved dielectric bodies such as spherical colloids or proteins of cylindrical geometry. The 
consideration of curved geometries deserve further investigation, since the curvature is known to come into play in 
various problems of biological interest, such as ionic exclusion from membrane nanopores of cylindrical geometry |44l 
145] or counterion condensation onto charged proteins 56]. Furthermore, the derived variational equations generalize 
the approximative self-consistent equations used in nanofiltration theories [50 - 52J to the case of ions associated with 
important core interactions and hydrodynamic forces deserve to be included in the model for applications related to 
ionic transport phenomenon |57j . 

The present model has of course several limitations. First of all, our variational calculation scheme is based on a 
generalized Onsager-Samaras approximation that assumes a uniform screening parameter for electrostatic and core 
interactions. As mentioned above, this simple choice might be responsible for the deviation of the variational result 
from simulation data in the close neighborhood of the interface. A more sophisticated variational choice, such as 
a spatially varying screening length is expected to correct the observed discrepancy. We currently investigate this 
point. Moreover, the variational equations were derived from a first order variational free energy, which limited our 
investigation to the submolar concentration regime. Corrections from a second order cumulant expansion might be 
necessary to consider higher concentration regimes. 

We finally note that in the case of ions confined into closed geometries characterized by a dielectric discontinuity, 
the existence of an infinite number of electrostatic images considerably complicates reliable numerical computations 
such as MC simulations. At the present, we are not aware of any simulation data for charged Yukawa liquids in 
contact with dielectric interfaces. The generalization of recent MC algorithms developed for confined Coulomb liquids 
to the case of charged Yukawa liquids will able us to determine in the future the validity domain of the present 
formalism EHl [Ml . 
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Appendix A: Boundary conditions for the external fields 



This Appendix is devoted to the derivation of the boundary conditions that come into play when solving Eqs. ( |31[ ) 
and (32) or Eqs. (31 1 and (32). In the case of a charged HC liquid confined between two planar surfaces located at 
z = and z = d, with surface charge a{z) = as [S{z) + 6{z — d)], the first step to derive the boundary conditions 
consists in integrating Eqs. (31) and (32 1 over the intervals z e [0^,0+] and z £ [a~,a+]. One obtains 







d(l) 




dz 


z=0+ 


dz 


2=0- 


d(t) 




d(j) 




dz 


z—a+ 


dz 


z—a^ 


dip 




dip 




dz 


z=0+ 


dz 


z=0~ 


dip 




dip 




dz 


2— a+ 


dz 


z—a^ 



Furthermore, in the particle free region z < a^, Pb = and Eqs. (31), (32) reduce to 





d^ 
dz^ 
d^iP 
'd^ 



b^ip = 0. 



The solution of Eqs. (A5) and (A6) that satisfy the boundary conditions < 

(p{z)=0, z<0 

(p{z) = Ci + C2Z , < z < 



-00) — and ip{z -t- —00) 



ip{z) = 036^^ , z < a. 



(Al) 
(A2) 
(A3) 
(A4) 

(A5) 
(A6) 

is 

(A7) 
(A8) 
(A9) 



where ci, 02,03 are constants of integ rati on. Deriving Eq. (A9), we assumed the continuity of the Yukawa field, i.e. 
■0(0"'") — ip{0~). According to Eq. (A7), the derivative of the electrostatic potential in the substrate part z < 
vanishes, d(p{z)/dz = 0, which yields with Eq. (Al) the Gauss law d(p{z) / dz]^^ — —Anqt^as- Thus Eq. (A8) becomes 

(p{z) = ci — AirqlwOsZ , < z < a^. (AlO) 



By injecting Eq. ( |A10| into Eq. ( [A2| ) and Eq. (|A9| into Eq. ( [Ai] ) and imp osing th e continuity of the Yukawa potential 
ip{z — a^) = ip{z = a^), we finally obtain the boundary conditions Eqs. (33)-(34) for the derivatives of the potentials 
at z = a„, . Following the same steps as above for the interface located at z = d — a^, , one gets the additional boundary 
conditions Eqs. (35)-(36). 



Appendix B: Variational Yukawa and electrostatic Green's functions 

In this Appendix, we explain the derivation of the variational HC potential wq{v^v') and the electrostatic potential 
uo(r,r') for the slit pore geometry depicted in Fig. [T] These potentials are obtained by inverting the relations (56) 
and (57), i.e. by solving a generalized DH equation 



[-V(e(r)V) + e(r)K2(r)] [/(r, r') = \5{t ~ r') 



(Bl) 



where the free parameter A will be adjusted later in order to recover wo(r,r') and K;o(r,r') from the kernel f7(r,r'). 
The piecewise dielectric permittivity and the inverse screening length, respectively, are given by 



£(z) = £>6i(z)6i(rf-z)+e<[6i(-z) + 6l(z-d)] 
Kj,(z) = K< [6'(a^ — z) + 0(z + flu, — d)] + K>6'(z — a,iu)0(d — fltu — z). 



By injecting into Eq. (Bl ) the Fourier expansion of the kernel 

l2l 



C/(r,r') = / 



47r2 



[/(z,z',ft). 



(B2) 
(B3) 



(B4) 
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one gets 



[-d,e{z)d, + e{z) {k^ + k^{z))] U{z, z' , k) = XS{z - z'). 



(B5) 



The homogeneous solution to this equation is of the form U = C+e^^^^''"^^ + C_e Vfe^+K^^^ ^pj-jg coefficients C± for 
each layer are found by imposing the finiteness of the kernel lim2_>.±oo U{z) — and the boundary conditions 



dU 



9z 



9z 
dz 



£(Z = I]+) 



5z 



(B6) 
(B7) 



(B8) 



where S denotes the position of the surfaces at z = 0^aw,d — aw,d and we also introduced the notation U {'E±) 
lime_^o f/ (S ± e). The kernel reads for < z,z' < d 



U{r,r') 



\ g-«>|r-r'| 

47re> |r — r'l 



SU{r,r') 



with the anisotropic part 
d^k 



SU{r,r') = 

5U{z,z',k) = 



47r2 
A 



k-(r-r') 



6U{z,z\k) 
A 



(B9) 



(BIO) 



2e>p> 1 - A2e-2p>(d-2a„) 



^p>(^+^'-2a„) _^gP>(z+.'+2a,„-2d) _^2Ae-2P>('i-2a™)cOSh (p>|z-z'|) 

(Bll) 



where we have defined 



A = , ri = 



P> + VP< ' 



1 + Aoe-2p<a„ 



(B12) 



and p< = ^/fc2 + K<, Aq = (e> — e<)/(e> + e<). By comparing the relations (561 and (57 1 with Eq. (Bll, the 



variational potentials wo(r,r') and wo(r,r') directly follow from the solution (Bll) as explained 



below. 



First, by setting k< = 0, k> = Kc, £< = £m: £> = £w, A = ineytyj and defining — \/ k'^ + k'^ 

Pc -Vck 1 - Aoe-2'^'^™ 



Pc + Vck ' 



1 + Aoe-2fea„ 



and Aq ~ {e^ — £m)/{£w +em), the variational electrostatic potential is obtained from Eq. (B9| in the form 

g-Kc|r-r'| 

vo{y,y') = £^,- — + (5wo(r,r') 

|r — r'l 



where 



<5t;o(r,r')= J ^e^'^<''-'''^ Svoiz, z' ,k) 



(B13) 



(B14) 



(B15) 



Svoiz, z', k) 



2-kL, 



1 - A2e-2Pc(d-2a„) 



g-p.(2+^'-2a„) ^gP,(^+2'+2a„-2d) _^2A,e-2p>('*-2'^™)cOsh(p>|z-z'|) . 

(B16) 



Second, by setting = = K,y, £< = e>, A — Ane^iy and defining pm — + fc^, Py — ^ i^y -\- kP- and 

^ Prn 



A„ 



(B17) 
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the Yukawa potential follows from Eq. ( B9 1 as 



g-K„|r-r'| 

wo{'r,r') =£y-^——^ + 5wo{r,r') 



where 



2-irf A 



(B18) 



(B19) 



1 _ A2p-2p„(d-2a„) 

y 



Py 



-Py(z+z' — 2a„) 



(B20) 



Appendix C: Computation of the surface tension with the charging procedure 



We explain in this Appendix the computation of the surface tension erg associated with a neutral interface. The 
surface tension is defined as the excess Grand potential, 



(Te = [fly - flvipb = 0) - rii] 

where the Grand potential of the bulk system given by Eq. ( [T5| reads 

[kI + {nyt - b){Kl^ + Kybb - 2b^)] - 



1 

24^ 



8ttI„ 



(CI) 



(C2) 



In the above equation, L stands for the length of the system. We also note that in Eq. (CI I, the subtracted quantity 
^o<p{Pb = 0) corresponds to the vdW level surface tension of the pure water. 

The derivation of the vdW part of the surface tension fio with the charging procedure was explained in Ref. [H] 
for the DH theory of slit pores. Using a general kernel as in Appendix [B] the Coulomb and Yukawa contributions 
Eqs. (46 1 and (47 1 will be derived with a single calculation. The charging method consists in reexpressing the excess 



van der Waals energy in the form of an integral over an auxiliary charging parameter ^, 



Aflo^ — r^oip - ^o^ipb = 0) = - In 



J Vip e 



(C3) 



de-ln / e-/ff.(r)[(V^)=+.i(r)] 





where the operator U ^ (r, r'; k{z)) is defined by Eq. (Bl ), whose solution for a single interface geometry follows from 



Eq. (BlOl in the limit d oo. In the same limit, the piecewise functions e{z) and k{z) defined in Eqs. (B2) and (B3) 
become = ey9{z) + e^9{—z) and Ky{z) — K^9{aw — z) + Ky9{z — ay,). As in Appendix pi the free parameter A will 
be adjusted later and the fluctuating poten tial (p will be identified with or ■!/; at the endof the calculation in order 
to recover ilo^ and flo^ from flo^. In Eq. (C3), we have also introduced the function K^(r) = + [^(r)^ — k^]. 
Evaluating the derivative with respect to ^ in Eq. ( C3 1 , one obtains 



AflQip — S 



£>(« 



2A 



d^ / dzU{r,r; K^) 



(C4) 



Carrying out the integral over z and substituting as in Appendix IB| A, Kg and eg by their values associated with the 



Yukawa and Coulomb contributions, one recovers Aflg^ and Af^o^ from Afig^. 



vdW-level contributions and the bulk Potential Eq. (C2 1 into the relation (CI I, one finally obtains the surface tension 
Eq. (iTl. 



Injecting ri^, of Eq. (581 with these 



Appendix D: Monte Carlo simulations 

We present in this Appendix the details of the canonical MC simulations. All simulation results presented in section 



IV B 



were obtained by performing standard Monte Carlo simulations |60) in a slab geometry of size Ix 



In 



-1/3 
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— 1/3 

and Iz = SOpj, , with impenetrable walls located at 2; = and z = Iz- Periodic boundary conditions in the x and 
y directions were used. The value of Iz that we considered is large enough so that the box geometry is equivalent to 
the single interface case presented in section [IV B| During the simulations, the total number of particles was fixed at 
Np = 20000, with the corresponding particle density pMC — Np/ [l^lylz) = Pb- The average acceptance rate during the 
simulations was approximately 0.80. For each set of parameters, the system is initialized with a random configuration 
and first equilibrated over 50000 Monte Carlo steps. After equilibration, the density profiles are obtained according 
to 

— 1/3 

where Zi is the coordinate of the particle i along z axis, dz ~ / 50 is the width of a transversal bin and 

< • > denotes the ensemble average. For all cases, the average was obtained from 20 independent runs with 5000 
measurements taken at 20 Monte Carlo steps interval for each run. The density profiles were obtained from an average 
over 20 independent runs, with 5000 measurements taken at intervals of 20 Monte Carlo steps for each run. 
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FIG. 2: (Color online) Density profile of neutral Yukawa particles at the interface for Xy = 2.992 and (a) Fy — 0.25, (b) 
Ty — 0.4 and (c) Ty = 0.6. The plots compare MC results (squares) with the variational method (continuous lines), the MF 
theory (dash-dotted line) and the RPA method (dashed lines). The black circles in (a) are obtained with the linear solution 
Eq. (1411. 
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FIG. 3: (Color online) Density profile of monovalent ions {q — 1) at the dielectric interface for a„ = 0, Eu, ~ 78, = 2 and 
ion sizes at = 0, 1.5,2,2.5,3 A. Bulk concentrations are (a) pb = 0.3 M and (b) pi, = 1 M. The dashed curves in the bottom 
plot shows the density profiles for = 1.5 A and = 3.0 A for the case Sm = £w = 78. 
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FIG. 4: (Color online) Solid lines illustrate for pt = 1 M, ai = 3.0 A and a„ = the behaviour of the potentials tp{z) — tpi,, 
Vy{z), Vcs{z), Vci{z) and the total PMF in units of fcsT. Brown and red dashed lines show 4}{z) — and Vy{z), respectively, 
for a smaller ion size of — 1.5 A. 
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FIG. 5: (Color online) Surface tension of monovalent ions against the bulk concentration for various ion sizes, Em = 2, (a) 
= and (b) a™ = Oi. 
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FIG. 7: (Color online) (a) Partition coefficient of neutral particles of size ai = 2 K against d for = at p(, = 0.2 M 
(red curves) and pb ~ 0.5 M (black curve). The dashed line is the prediction of the MF theory and the solid lines display the 
variational result, (b) Variational inverse screening length Hy vs d for neutral particles (blue dotted-dashed line) and monovalent 
ions in the cases Em = (solid black line) and Sm = 2 (solid red line), pb = 0.5 M and the other model parameters are the 
same as in (a). 
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(b) 

FIG. 8: (Color online) Density profile of neutral Yukawa particles (blue dash-dotted line), monovalent ions without size (dashed 
lines) and with finite size = 2 A (solid lines) for (a) d = 8 A and (b) d = 16 A. Model parameters = A, p6 = 0.5 

M, Em = 78 (black lines) and Sm = 2 (red lines). 
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FIG. 9: (Color online) (a) Partition coefficient of monovalent ions without size (dashed lines) and of finite size = 2 A (solid 
lines) against the pore size d for a^j = A, pb — 0.5 M, e„ = 78 and £m = 2 to 78. (b) Partition coefficient against pb with 
d — 10 K and the same model parameters as in (a). 
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FIG. 10: (Color online) Variational inverse screening length acc against the pore size d for monovalent ions without size (dashed 
lines) and with finite size ai — 2 A (solid lines). Model parameters are = A, pb — 0.5 M, Ew = 78 (black lines) and Em = 2 
(dashed lines). 
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FIG. 11: (Color online) Partition coefficient of monovalent ions against the pore size d with a finite Stern layer of thickness 
— Ui. Model parameters are the same as in Fig. |9]a. 



